---
title: "Replication of Results of 'Social Trust and the Winner-Looser Gap'"
author: "Matias Bargsted \& Andrés González-Ide"
date: "2023-06-09"
output:
  pdf_document:
    number_section: true
header-includes:
    - \usepackage{hyperref} 
    - \usepackage{float}
    - \floatplacement{figure}{H}
    - \renewcommand{\thesection}{\Alph{section}}
    - \usepackage{caption}
    - \DeclareCaptionLabelFormat{nospace}{#1#2}
    - \captionsetup[figure]{labelformat=nospace}
    - \captionsetup[table]{labelformat=nospace}
bibliography: "https://api.citedrive.com/bib/41a0353d-291a-49e7-b197-43f2f2e02443/references.bib?x=eyJpZCI6ICI0MWEwMzUzZC0yOTFhLTQ5ZTctYjE5Ny00M2YyZjJlMDI0NDMiLCAidXNlciI6ICI0ODk3IiwgInNpZ25hdHVyZSI6ICI0MmZkYzFlZmU3YTE1ZjAwNTQ3OWQ5NzZlNTkwYjNkM2NkNzY4NjU1NmQ4ZWU5YjA1MmY4NDE2N2U1MTE2YmFhIn0=/bibliography.bib"
---

```{r database work 2021, include = FALSE}

knitr::opts_chunk$set(echo = FALSE,warning = FALSE,message = FALSE,results = "asis")
library(stargazer)
library(plm)
library(texreg)
library(AER)
library(DRDID)
library(survey)
library(questionr)
library(sjmisc)
library(dplyr)
#library(seres)
library(gurobi)
library(MatchIt)
library(cobalt)
library(gridExtra)
library(ggplot2)
library(svrep)
library(haven)
library(Hmisc)
library(xtable)
library(psych)
library(lavaan)
library(semPlot)
library(ggcorrplot)

ruta <- "C:/Users/Matias/Dropbox/Proyectos/Fondecyt Participacion Electoral/Articulos/Electoral Winners y Plebiscito/Harvard Dataverse/"
#ruta <- "C:/Users/aigon/Dropbox/Fondecyt Participacion Electoral/"

#casen <- read_dta(paste0(ruta,"Encuesta Plebiscito 2022/3_Datos/CASEN 2020/CASEN 2020.dta"))
#casen <- casen %>% filter(edad %in% c(18:65))

sum.func <- function(x, w) {
  # Verify dichotomous values
  unique_vals <- unique(na.omit(x))
  is_dichotomous <- length(unique_vals) == 2 && all(unique_vals %in% c(0, 1, 2))
  
  # Summary statistics
  if (is_dichotomous) {
    cbind(table(is.na(x))[1],
          round(mean(x, na.rm=TRUE), 2),
          "-", # SD
          round(weighted.mean(x, w, na.rm=TRUE), 2),
          "-", # WSD
          round(min(x, na.rm=TRUE), 2),
          round(quantile(x, probs=0.25, na.rm=TRUE), 2),
          round(quantile(x, probs=0.75, na.rm=TRUE), 2),
          round(max(x, na.rm=TRUE), 2))
  } else {
    cbind(table(is.na(x))[1],
          round(mean(x, na.rm=TRUE), 2),
          round(sd(x, na.rm=TRUE), 2),
          round(weighted.mean(x, w, na.rm=TRUE), 2),
          round(sqrt(Hmisc::wtd.var(x, w, na.rm=TRUE)), 2),
          round(min(x, na.rm=TRUE), 2),
          round(quantile(x, probs=0.25, na.rm=TRUE), 2),
          round(quantile(x, probs=0.75, na.rm=TRUE), 2),
          round(max(x, na.rm=TRUE), 2))
  }
}


#sum_vars <- with(pb, data.frame(
#  confianza_w01 = conf_soc_w01, confianza_w02 = conf_soc_w02,
#  sex = (as.numeric(mujer)-1), age = edad, education = as.numeric(educa_w01), 
#  mean_years_education = anio_educ_prom_w01
#  ))

#sum_vars <- na.omit(sum_vars)

```


```{r database work,results='hide'}

###PRESIDENTIAL ELECTION STUDY

#Open files
load(paste0(ruta,"1. Presidential Election (study 1)/2. Database/MAYO_WIDE_RECODE.Rdata"))
mayo <- mayo_wide; mayo_wide <- NULL

## Procesamiento de datos

table(as.numeric(mayo$votacion_segunda_w04))
table(mayo$voto_2vuelta_w05)
mayo$delta_elec_winner <- as.numeric(mayo$voto_2vuelta_w05=="Gabriel Boric")-0
mayo$delta_elec_winner <- ifelse(mayo$voto_2vuelta_w05 %in% c("Jos? Antonio Kast","Gabriel Boric"),mayo$delta_elec_winner,NA)
table(mayo$delta_elec_winner)

#######---------- Creacion de Variables --------------##########

#Edad c("18 a 29","30 a 39","40 a 49", "50 a 65")
mayo$edad_w04F <- factor(car::recode(mayo$edad_w04, "18:29=1; 30:39=2; 40:49=3; 50:65=4"),
                           labels=c("18-29","30-39","40-49", "50-65"))
prop.table(table(mayo$edad_w04F))

#Sexo
mayo$mujer_w03 <- as.numeric(mayo$sexo_w03==2)
prop.table(table(mayo$mujer_w03))

#Educacion
mayo$educaF_w03 <- factor(car::recode(mayo$educ_w03, "1:5=1; 6:7=2; 8:10=3; 99=1"),
                            labels=c("Secondary Ed.", "Tertiary Tech.", "Tertiary Univ"))
prop.table(table(mayo$educaF_w03))

#Political Id. 

mayo$lrscaleF_w04 <- factor(car::recode(mayo$P15_SQ001_w04, 
                                          "0:4=1; 5=2; 6:10=3;else=4"),
                              labels = c("Izquierda", "Centro", "Derecha", "No Identificado"))
table(mayo$lrscaleF_w04)

#Interest in politics
mayo$pol_interesF_w04 <- as.numeric(mayo$pol_interes_w04)>=4
prop.table(table(mayo$pol_interesF_w04))
mayo$pol_interesF2_w04 <- factor(case_when(
  as.numeric(mayo$pol_interes_w04) %in% c(1,2) ~ 1,
  as.numeric(mayo$pol_interes_w04) == 3 ~ 2,
  as.numeric(mayo$pol_interes_w04) %in% c(4,5) ~ 3
),labels = c("Sin interes","Algo interesado","Interesado"))

#Create Dataset with key variables only and no missing (to help id outliers)
keep <- c("conf_soc_w05", "conf_soc_w04", "conf_soc_w03", "delta_elec_winner", "pol_interesF_w04",
          "pol_interesF2_w04",
          "lrscaleF_w04","educ_w03","educ_w01", "educaF_w03", "edad_w01", "edad_w04F", "mujer_w03",
          "atricion_w05","atricion_w04",
          "conf1_w05","conf2_w05","conf3_w05","conf1_w04","conf2_w04","conf3_w04",
          "conf1_w03","conf2_w03","conf3_w03")
mm <- na.omit(subset(mayo, select=keep))
rownames(mm) <- 1:nrow(mm)

alpha(with(mm,cbind(conf1_w05,conf2_w05,conf3_w05)))
alpha(with(mm,cbind(conf1_w04,conf2_w04,conf3_w04)))


```


```{r database work 2022, include=FALSE}

#### Plebiscito 2022

load(paste0(ruta,"2. Constitutional Plebiscite (study 2)/2. Database/Encuesta PB-2022.Rdata"))
ls()
names(pb2022)
pb2022 <- pb2022 %>% filter(attribute_235_w01 %in% c(18:65))


#######---------- Creacion de Variables --------------##########

#Electoral winning status
##Substract treatment in time 2 minus treatment in time 1 (when no one is treated)
table(pb2022$voto_sept_w02)
pb2022$delta_elec_winner <- as.numeric(pb2022$voto_sept_w02=="Rechazo") - 0 
prop.table(table(pb2022$delta_elec_winner))
pb2022$delta_elec_winner <- ifelse(pb2022$voto_sept_w02=="Nulo/Blanco",NA,pb2022$delta_elec_winner)
prop.table(table(pb2022$delta_elec_winner))

#Edad
pb2022$edad_w01 <- pb2022$attribute_235_w01
pb2022$edad_w01F <- factor(car::recode(pb2022$edad_w01, "18:29=1; 30:39=2; 40:49=3; 50:65=4"),
                           labels=c("18-29","30-39","40-49", "50-65"))
prop.table(table(pb2022$edad_w01F))

#Sexo
pb2022$mujer_w01 <- as.numeric(pb2022$attribute_224_w01==2)
prop.table(table(pb2022$mujer_w01))

#Educacion
pb2022$educaF_w01 <- factor(car::recode(pb2022$educ_w01, "1:5=1; 6:7=2; 8:10=3; 99=1"),
                            labels=c("Secondary Ed.", "Tertiary Tech.", "Tertiary Univ"))
prop.table(table(pb2022$educaF_w01))

#Diference Const & Person (ideological poostion)
pb2022$dist_const_pers_w01 <- (pb2022$lrscaleC_cont_w01 - pb2022$lrscaleC_w01)^2 

#Tramos expectativas
cutpoints = seq(min(pb2022$expect_apruebo_w01),max(pb2022$expect_apruebo_w01),5)
pb2022$expect_apruebo_w01R = NA
pb2022$expect_apruebo_w01R[pb2022$expect_apruebo_w01<cutpoints[1]] = 1
for(i in 1:length(cutpoints)) {
  pb2022$expect_apruebo_w01R = ifelse(pb2022$expect_apruebo_w01>=cutpoints[i] & pb2022$expect_apruebo_w01<cutpoints[i+1], i+1, 
                                      pb2022$expect_apruebo_w01R)
}

pb2022$expect_rechazo_w01 = NA
pb2022$expect_rechazo_w01 = 100-pb2022$expect_apruebo_w01

cutpoints = seq(min(pb2022$expect_rechazo_w01),max(pb2022$expect_rechazo_w01),5)
pb2022$expect_rechazo_w01R = NA
pb2022$expect_rechazo_w01R[pb2022$expect_rechazo_w01<cutpoints[1]] = 1
for(i in 1:length(cutpoints)) {
  pb2022$expect_rechazo_w01R = ifelse(pb2022$expect_rechazo_w01>=cutpoints[i] & pb2022$expect_rechazo_w01<cutpoints[i+1], i+1, 
                                      pb2022$expect_rechazo_w01R)
}

pb2022$expect_apruebo_w01R <- factor(case_when(
  pb2022$expect_apruebo_w01R %in% c(1:11) ~ 1,
  pb2022$expect_apruebo_w01R %in% c(12:15) ~ 2,
  pb2022$expect_apruebo_w01R %in% c(16:21) ~ 3),labels = c("5-50","55-70","75-100"))

pb2022$expect_rechazo_w01R <- factor(case_when(
  pb2022$expect_rechazo_w01R %in% c(1:11) ~ 1,
  pb2022$expect_rechazo_w01R %in% c(12:15) ~ 2,
  pb2022$expect_rechazo_w01R %in% c(16:21) ~ 3),labels = c("5-50","55-70","75-100"))

pb2022$dif_apbrech = NA
pb2022$dif_apbrech = abs(pb2022$expect_rechazo_w01-pb2022$expect_apruebo_w01)
cutpoints = seq(min(pb2022$dif_apbrech),max(pb2022$dif_apbrech),5)
pb2022$dif_apbrechR = NA
pb2022$dif_apbrechR[pb2022$dif_apbrech<cutpoints[1]] = 1
for(i in 1:length(cutpoints)) {
  pb2022$dif_apbrechR = ifelse(pb2022$dif_apbrech>=cutpoints[i] & pb2022$dif_apbrech<cutpoints[i+1], i+1, 
                                      pb2022$dif_apbrechR)
}

#Left-right scale grouped
pb2022$lrscaleF_w01 <- factor(car::recode(pb2022$lrscaleC_w01, 
                                          "0:4=1; 5=2; 6:10=3;else=4"),
                              labels = c("Izquierda", "Centro", "Derecha", "No Identificado"))
table(pb2022$lrscaleF_w01)
table(pb2022$delta_elec_winner, pb2022$lrscaleF_w01)

#A previous favorable opinion towards the act of voting is expected to magnify
#the effect of preference agreement between voters and the electorate. Individuals who 
#believe that the act of voting influences the outcome will see corroborated their prior 
#belief, which would 

#Interest in politics
pb2022$pol_interesF_w01 <- as.numeric(pb2022$pol_interes_w01)>=4
pb2022$pol_interesF2_w01 <- factor(case_when(
  as.numeric(pb2022$pol_interes_w01) %in% c(1,2)~1,
  as.numeric(pb2022$pol_interes_w01) == 3 ~ 2,
  as.numeric(pb2022$pol_interes_w01) %in% c(4,5)~3
),labels=c("Sin interes","Algo de interes","Interesado"))

prop.table(table(pb2022$pol_interesF_w01))

#Education as predictor of Engagement & Social Trust
prop.table(table(pb2022$pol_interesT_w01, pb2022$educaF_w01),2)

```


```{r pre-match weightning presidential,results='hide',fig.show='hide'}
## PRE-MATCH WEIGHTNING - PRESIDENTIAL ELECTION SURVEY
source(paste0(ruta,"3. Analysis Code/Weight_Presidencial_Pre-Match_v2.R"))
summary(mm$weight.educ.edad.pre_match)

```

```{r pre-match weighting 2022, results='hide',fig.show='hide'}
## PRE-MATCH WEIGHTNING -  PLEBISCITE SURVEY
source(paste0(ruta,"3. Analysis Code/Weight_Plebiscito_Pre-Match_v2.R"))
summary(pb2022$weight.educ.edad.pre_match)

```


```{r alphas, results='asis'}

# Presidential election

alfa_mm1 <- alpha(with(mm,cbind(conf1_w04,conf2_w04,conf3_w04)))
alfa_mm2 <- alpha(with(mm,cbind(conf1_w05,conf2_w05,conf3_w05)))

# Constitutional plebiscite

alfa_pb1 <- alpha(with(pb2022,cbind(P18_1_w01,P18_2_w01,P18_3_w01)))
alfa_pb2 <- alpha(with(pb2022,cbind(P10_1_w02,P10_2_w02,P10_3_w02)))

alphas <- rbind(alfa_mm1$total,alfa_mm2$total,alfa_pb1$total,alfa_pb2$total)
row.names(alphas) <- c("Pre-Presidential \nElection","Post-Presidential \nElection",
                       "Pre-Constitutional \nPlebiscite","Post-Constitutional \nPlebisicte")

colnames(alphas) <- c("Raw Alphas","Std.Aplha","G6 (smc)",
                      "Average R","S/N","ASE","Mean","SD","Median")

print(xtable(alphas,align=c("l",rep("c",9)),caption = "Cronbachs' Alphas for Social Trust for the Presidential and Constitutional Plebisicte Surveys",
             label = "table:alphas"),scalebox = .85,table.placement = "H",comment = F)

```


# Summary Statistics.

\setcounter{table}{0}
\def\figurename{Figure B}
\def\tablename{Table B}

The Tables B\ref{table:2021} and B\ref{table:2022} depict the summary statistics of all the variables employed in the study. The Population column shows the proportion of the population of each socio-demographic variable used in the study according to the 2020 Socio-economic Characterization Survey (CASEN). The CASEN survey is a large household probabilistic survey carried out by the Chilean Ministry of Social Development and Family periodically since 1990. The sample of the 2020 wave contained 185,437 respondents. The main objective of this survey is to measure the distribution of income among the population, but it also contains a rich battery of socio-demographic variables, such as education, marital status, household composition, well-being, among others. The full data set and information can be found at \url{https://observatorio.ministeriodesarrollosocial.gob.cl/encuesta-casen}. The Population column also shows the official results for each election according to the Chilean Electoral Service (SERVEL).

```{r summary table, results = 'asis'}

# CASEN 2020 Data

#casen2020 <- casen[c("educ","edad","sexo","expr")]
#casen2020 <- casen2020 %>% filter(edad %in% c(18:65), educ != 99) %>% mutate(
#  edad_r = factor(case_when(
#    edad %in% c(18:29) ~ 1,
#    edad %in% c(30:39) ~ 2,
#    edad %in% c(40:49) ~ 3,
#    edad %in% c(50:65) ~ 4
#  ),labels = c("18-29","30-39","40-49","50-65")),
#  educ_r = factor(case_when(
#    educ %in% c(0:6) ~ 1,
#    educ %in% c(7,8) ~ 2,
#    educ %in% c(9:12) ~ 3
#  ), labels = c("Primary and Secondary Education",
#  "Tertiary Technical Education","Tertiary Professional Education")),
#  educ_r2 = factor(case_when(
#    educ %in% c(0:6) ~ 1,
#    educ == 7 ~ 2,
#    educ == 8 ~ 3,
#    educ == 9 ~ 4,
#    educ %in% c(10:12) ~ 5
#  ), labels = c("Primary and Secondary Education",
#  "Tertiary Technical Education Incomplete","Tertiary Technical Education Complete","Tertiary Professional Education Incomplete","Tertiary Professional Education Complete")))

## 2021

sum_vars2021 <- with(mm,data.frame(voto_boric = delta_elec_winner,
  social_trust_w04 = conf_soc_w04, social_trust_w05 = conf_soc_w05,
  sex = as.numeric(mujer_w03), age = edad_w04F, education = educaF_w03, 
  pol_interes = as.numeric(pol_interesF2_w04),izq_der = lrscaleF_w04))

sum_vars2021 <- na.omit(sum_vars2021)

sum_vars2021$Center <- as.numeric(sum_vars2021$izq_der == "Centro")
sum_vars2021$Right <- as.numeric(sum_vars2021$izq_der == "Derecha")
sum_vars2021$Left <- as.numeric(sum_vars2021$izq_der == "Izquierda")
sum_vars2021$NoId <- as.numeric(sum_vars2021$izq_der == "No Identificado")
sum_vars2021$izq_der <- NULL

sum_vars2021$low_pi <- as.numeric(sum_vars2021$pol_interes==1)
sum_vars2021$mid_pi <- as.numeric(sum_vars2021$pol_interes==2)
sum_vars2021$hig_pi <- as.numeric(sum_vars2021$pol_interes==3)
sum_vars2021$pol_interes <- NULL

sum_vars2021$age18 <- as.numeric(sum_vars2021$age == "18-29")
sum_vars2021$age30 <- as.numeric(sum_vars2021$age == "30-39")
sum_vars2021$age40 <- as.numeric(sum_vars2021$age == "40-49")
sum_vars2021$age50 <- as.numeric(sum_vars2021$age == "50-65")
sum_vars2021$age <- NULL

sum_vars2021$secondary <- as.numeric(sum_vars2021$education == "Secondary Ed.")
sum_vars2021$tertiary.t <- as.numeric(sum_vars2021$education == "Tertiary Tech.")
sum_vars2021$tertiary.u <- as.numeric(sum_vars2021$education == "Tertiary Univ")
sum_vars2021$education <- NULL

sum_vars2021 <- sum_vars2021 %>% relocate(sex,.before = age18)

sum.data21 <- data.frame(t(sapply(1:dim(sum_vars2021)[2],
                                  function(x) sum.func(sum_vars2021[mm$atricion_w05=="Permanece",x],
                                                       mm$weight.educ.edad.pre_match[mm$atricion_w05=="Permanece"]))))

sum.data21$casen <- "-"

rownames(sum.data21) <- c("Vote for Boric","Social Trust W04","Social Trust W05",
                          "Center","Right","Left","No Political Id.","Low Political Interest",
                          "Medium Political Interest","High Political Interest","Female",
                          "18-29","30-39","40-49","50-65","Secondary Ed.","Tertiary Tech.",
                          "Tertiary Univ.")


colnames(sum.data21) <- c("N", "Mean", "St. Dev.", "Weighted Mean", "Weighted St. Dev.",
                        "Min", "Pctl(25)", "Pctl(75)", "Max","Population")

#casen <- casen %>% mutate(edad_r = factor(case_when(
#  edad %in% c(18:29) ~ 1,
#  edad %in% c(30:39) ~ 2,
#  edad %in% c(40:49) ~ 3,
#  edad %in% c(50:65) ~ 4
#),labels = c("18-29","30-39","40-49","50-65")))
#prop.table(questionr::wtd.table(casen$edad_r,weights = casen$expr))*100

sum.data21$Population[which(row.names(sum.data21)=="Vote for Boric")] <- 0.56
sum.data21$Population[which(row.names(sum.data21)=="Female")] <- 0.55
sum.data21$Population[which(row.names(sum.data21)=="18-29")] <- 0.29
sum.data21$Population[which(row.names(sum.data21)=="30-39")] <- 0.21
sum.data21$Population[which(row.names(sum.data21)=="40-49")] <- 0.19
sum.data21$Population[which(row.names(sum.data21)=="50-65")] <- 0.31
sum.data21$Population[which(row.names(sum.data21)=="Secondary Ed.")] <- 0.58
sum.data21$Population[which(row.names(sum.data21)=="Tertiary Tech.")] <- 0.13
sum.data21$Population[which(row.names(sum.data21)=="Tertiary Univ.")] <- 0.29

## 2022

sum_vars2022 <- with(pb2022, data.frame(voto_rechazo = delta_elec_winner,
  social_trust_w01=prs_trst_w01, social_trust_w02 = prs_trst_w02,
  sex = as.numeric(attribute_224_w01) - 1, age = edad_w01F, education = educaF_w01, 
  izq_der = lrscaleF_w01, pol_interes = as.numeric(pol_interesF2_w01)
))

sum_vars2022 <- na.omit(sum_vars2022)

sum_vars2022$Center <- as.numeric(sum_vars2022$izq_der == "Centro")
sum_vars2022$Right <- as.numeric(sum_vars2022$izq_der == "Derecha")
sum_vars2022$Left <- as.numeric(sum_vars2022$izq_der == "Izquierda")
sum_vars2022$NoId <- as.numeric(sum_vars2022$izq_der == "No Identificado")
sum_vars2022$izq_der <- NULL

sum_vars2022$low_pi <- as.numeric(sum_vars2022$pol_interes==1)
sum_vars2022$mid_pi <- as.numeric(sum_vars2022$pol_interes==2)
sum_vars2022$hig_pi <- as.numeric(sum_vars2022$pol_interes==3)
sum_vars2022$pol_interes <- NULL

sum_vars2022$age18 <- as.numeric(sum_vars2022$age == "18-29")
sum_vars2022$age30 <- as.numeric(sum_vars2022$age == "30-39")
sum_vars2022$age40 <- as.numeric(sum_vars2022$age == "40-49")
sum_vars2022$age50 <- as.numeric(sum_vars2022$age == "50-65")
sum_vars2022$age <- NULL

sum_vars2022$secondary <- as.numeric(sum_vars2022$education == "Secondary Ed.")
sum_vars2022$tertiary.t <- as.numeric(sum_vars2022$education == "Tertiary Tech.")
sum_vars2022$tertiary.u <- as.numeric(sum_vars2022$education == "Tertiary Univ")
sum_vars2022$education <- NULL

sum_vars2022 <- sum_vars2022 %>% relocate(sex,.before = age18)

sum.data22 <- data.frame(t(sapply(1:dim(sum_vars2022)[2],
                                  function(x) sum.func(sum_vars2022[pb2022$atricion=="Permanece",x],
                                                       pb2022$weight.educ.edad.pre_match[pb2022$atricion=="Permanece"]))))

sum.data22$casen <- "-"

rownames(sum.data22) <- c("Vote Reject","Social Trust W04","Social Trust W05",
                          "Center","Right","Left","No Political Id.","Low Political Interest",
                          "Medium Political Interest","High Political Interest","Female",
                          "18-29","30-39","40-49","50-65","Secondary Ed.","Tertiary Tech.",
                          "Tertiary Univ.")

colnames(sum.data22) <- c("N", "Mean", "St. Dev.", "Weighted Mean", "Weighted St. Dev.",
                        "Min", "Pctl(25)", "Pctl(75)", "Max","Population")

sum.data22$Population[which(row.names(sum.data22)=="Vote Reject")] <- 0.62
sum.data22$Population[which(row.names(sum.data22)=="Female")] <- 0.55
sum.data22$Population[which(row.names(sum.data22)=="18-29")] <- 0.29
sum.data22$Population[which(row.names(sum.data22)=="30-39")] <- 0.21
sum.data22$Population[which(row.names(sum.data22)=="40-49")] <- 0.19
sum.data22$Population[which(row.names(sum.data22)=="50-65")] <- 0.31
sum.data22$Population[which(row.names(sum.data22)=="Secondary Ed.")] <- 0.58 
sum.data22$Population[which(row.names(sum.data22)=="Tertiary Tech.")] <- 0.13
sum.data22$Population[which(row.names(sum.data22)=="Tertiary Univ.")] <- 0.29 

print(xtable(sum.data21, digits=c(0, 0, rep(2, length(colnames(sum.data21))-1)),
             align=c("l","c","c","c", "p{1.5cm}","p{1.5cm}","c","c","c","c","c"),
             caption = "Summary Statistics for the 2021 Presidential Election Variables",
             label = "table:sum2021"),
      scalebox=.85, table.placement="H",comment = FALSE)

print(xtable(sum.data22, digits=c(0, 0, rep(2, length(colnames(sum.data22))-1)),
             align=c("l","c","c","c", "p{1.5cm}","p{1.5cm}","c","c","c","c","c"),
             caption = "Summary Statistics for the 2022 plebiscite Variables",
             label = "table:sum2022"),
      scalebox=.85, table.placement="H",comment = FALSE)
```

# Survey post-estratification weighting

\setcounter{table}{0}
\def\figurename{Figure C}
\def\tablename{Table C}

Post stratification weights were calculated using the 2020 CASEN survey as our source of population values. Given that our samples under-represent, relative to the population, younger and older people, and over-represent those with high educational attainment our weights were based on age and education. The weights were calculated using the \textit{rake} function from the \textit{survey} package version 4.1.1 [@LumleyMisc]. 

Previous to the calculation of weights, education and age contained within the CASEN survey were recoded as follows:

\begin{itemize}
\item Education: respondents education level was grouped into five categories: (1) Primary and Secondary Education; (2) Tertiary Technical Education Incomplete; (3) Tertiary Technical Education Complete; (4)  Tertiary University Education Incomplete; and (5) Tertiary University Education Complete.
\item Age: we grouped the respondents in four categories: (1) 18-29; (2) 30-39; (3) 40-49; and (4) 50-65.  
\end{itemize}

Table C\ref{table:casen} depicts the distribution of the categories among the population of both variables.

```{r casen table, results='asis'}
table.df <- data.frame(
  N = c(
    3553186.00, 2612533.00, 2395648.00, 3799486.00, 
    7155934.00, 547296.00, 1120183.00, 1347058.00, 2190382.00
  ),
  Freq = c(
    28.75, 21.14, 19.38, 30.74, 
    57.89, 4.43, 9.06, 10.90, 17.72
  )
)

row.names(table.df) <- c(
    "18-29", "30-39", "40-49", "50-65", 
    "Primary and Secondary Education", 
    "Tertiary Technical Education Incomplete", 
    "Tertiary Technical Education Complete", 
    "Tertiary Professional Education Incomplete", 
    "Tertiary Professional Education Complete"
  )

print(xtable(table.df, align = c("l","c","c"),caption = "Frquency Table of Age and Education in the Population",
             label = "table:casen"),scaleboox = .8,table.placement = "H",comment = F)
```

Survey weights were calculated so the sample distribution of education and age resembled the distribution captured by the CASEN 2020 survey. However, this led to some unusually large values ($w_i>6$) that had a deleterious effects on the sampling variances of our regression estimates. Following Potter & Zheng [-@potter2015methods] and Van de Kerckhove, Mohadjer & Krenzke [-@vandeKerckhove2014], we trimmed the weights using the Interquartile range procedure which entails constraining the values of the weights within a range with a minimum value of 1/R and a maximum value of R, where R is defined as:

\begin{equation}
R = Me(w_i) + k*IQR(w_i)
\end{equation}

and \textit{Me} corresponds to the Median of the weight, $k$ is a constant that Potter & Zheng [-@potter2015methods] recommend setting to 4 or 5, and IQR is the interquartile range of the weight. According to the authors, this trimming procedure is effective in reducing the sample variances inflated by extreme values of the weights, but they do not prove a rationale for setting the constant $k$ to 4 or 5. Given that this seems a bit arbitrary, we calculated the mean square error (MSE) of the weighted sample proportions of people with secondary and university level education at different values of $k$, and selected the value where further increases of $k$ implied an seemingly irrelevant reduction in the MSE. The results of this analysis are shown in Figures C1 (presidential election survey) and C2 (constitutional plebiscite survey). As can be seen in graph A of both figures the MSE decreases following a quadratic pattern, where further increases beyond $k=9$ reduce the MSE in an increasingly smaller magnitude, especially for the proportion of respondents with secondary education or less. This can be seen more explicitly in graph B of both figures which plots on the Y\-axis the change in MSE between the iterative value of $k$ and $k-1$. Based on these insights we choose to set $k=9$, although slightly lower values could also be reasonable options. 
 

```{r 2021 mse plot,results='asis',fig.cap="\\label{fig:mse_pres}MSE reduction by K for the 2021 Presidential Election",fig.pos="H"}
gridExtra::grid.arrange(mse_plot21,delta_plot21)
```


```{r 2022 mse plot, results='asis',fig.cap="\\label{fig:mse_2022}MSE reduction by K for the 2022 Constitutional Plebiscite",fig.pos="H"}
gridExtra::grid.arrange(mse_plot22,delta_plot22)
```


The following tables depicts the distribution of the weights in both samples:

```{r summary weight pre-match, results='asis'}
stargazer(mm[c("weight.educ.edad.pre_match")],summary=T,header=F,
          title = "Summary statistics of 2021 presidential election survey weights",
          table.placement = "h",
          covariate.labels = c("Educational and Age Weight"))

stargazer(pb2022[c("weight.educ.edad.pre_match")],summary=T,header=F,
          title = "Summary statistics of 2021 constitutional plebiscite survey weights",
          table.placement = "h",
          covariate.labels = c("Educational and Age Weight"))

```

Equivalent weights were constructed for the subsample of matched cases of both studies. These were calculated using the same population proportions as reference and value of $k$ for weight trimming.

\newpage

# Matching pre-processing.

\setcounter{table}{0}
\setcounter{figure}{0}
\def\figurename{Figure D}
\def\tablename{Table D}

The sample of respondents of both surveys were matched using cardinality matching based on the same set of pre-treatment (pre-election wave) covariates: (a) ideological position (recoded in 4 categories), (b) Educational attainment (recoded in 3 categories), (c) Age (recoded in 3 categories), (d) Sex (binary), and (e) Political interest (recoded in 3 categories). The imbalance tolerance for mean differences for all the covariates in both surveys was set to 0.1. The solver used was Gurobi from the \textit{gurobi} R package version 10.0.1 [@gurobi]. The matched samples contain 380 cases for the 2021 presidential election survey, and 418 cases for the 2022 constitutional plebiscite survey. Figures D\ref{fig:2021} and D\ref{fig:2022} illustrate the differences among electoral winners and losers in the matched and pre-matched samples for the presidential election and for the constitutional plebiscite, respectively. As can be seen, the matching processing reduces drastically the pre-treatment differences among the groups in both samples, specially in what respects to ideological self-positioning and age. 

```{r matching 2021,results='asis',fig.cap="\\label{fig:2021}Matched and Unmatched Samples for 2021 Presidential Election Survey"}
#Matching Corruption Treatment 

m1.full <- matchit(delta_elec_winner ~ lrscaleF_w04 + educaF_w03 + edad_w04F + mujer_w03 + 
                     pol_interesF2_w04, data=mm, 
                   method = "cardinality", tols = .1, solver="gurobi")
#summary(m1.full)
#with(match.data(m1.full), table(delta_elec_winner,lrscaleF_w04))

## Cambiar los nombres

v <- data.frame(old = c("lrscaleF_w04_Izquierda","lrscaleF_w04_Centro","lrscaleF_w04_Derecha",
                        "lrscaleF_w04_No Identificado","educaF_w03_Secondary Ed.",
                        "educaF_w03_Tertiary Tech.",
                        "educaF_w03_Tertiary Univ","edad_w04F_18-29","edad_w04F_30-39","edad_w04F_40-49","edad_w04F_50-65",
                        "mujer_w03","pol_interesF2_w04_Sin interes","pol_interesF2_w04_Algo interesado","pol_interesF2_w04_Interesado",
                        "distance"),
                new = c("Left", "Center", "Right", "No Political Id.", 
                        "Secondary Ed.", "Tertiary Tech.", "Tertiary Univ.", "18-29","30-39","40-49","50-65",
                        "Women","Low Political Interest","Medium Political Interest","High Political Interest","Propensity Score"))


#summary(m1.full)

love.plot(m1.full, abs=T, title="Covariate Balance for Electoral Winning (2021)",
                        sample.names=c("Unmatched", "Matched"), stars = "std",
                        var.names = v) + theme_bw() + theme(legend.position = "bottom")

mayo.mt <- match.data(m1.full)
row.names(mayo.mt) <- 1:nrow(mayo.mt)

```


```{r matching 2022,results='asis',fig.cap="\\label{fig:2022}Matched and Unmatched Samples for 2022 Constitutional Plebiscite Survey",fig.pos="h"}
#### 2022 ####

m1.full.22 <- matchit(delta_elec_winner ~ lrscaleF_w01 + educaF_w01 + edad_w01F + mujer_w01 + 
                     pol_interesF2_w01, data=subset(pb2022, !is.na(delta_elec_winner)), 
                   #method = "full", replace=T)
                   method = "cardinality", tols = .1, solver="gurobi")

v22 <- data.frame(old = c("lrscaleF_w01_Izquierda","lrscaleF_w01_Centro","lrscaleF_w01_Derecha",
                        "lrscaleF_w01_No Identificado","educaF_w01_Secondary Ed.",
                        "educaF_w01_Tertiary Tech.",
                        "educaF_w01_Tertiary Univ","edad_w01F_18-29","edad_w01F_30-39","edad_w01F_40-49","edad_w01F_50-65",
                        "mujer_w01","pol_interesF2_w01_Sin interes","pol_interesF2_w01_Algo de interes","pol_interesF2_w01_Interesado",
                        "distance"),
                new = c("Left", "Center", "Right", "No Political Id.", 
                        "Secondary Ed.", "Tertiary Tech.", "Tertiary Univ.", "18-29","30-39","40-49","50-65",
                        "Women","Low Political Interest",
                        "Medium Political Interest","High Political Interest","Propensity Score"))
love.plot(m1.full.22, abs = T, title = "Covariate Balance for Electoral Winning (2022)",
                     sample.names = c("Unmatched","Matched"),stars="std",
                     var.names=v22) + theme_bw() + theme(legend.position = "bottom")

pb2022mt <- match.data(m1.full.22)
row.names(pb2022mt) <- 1:nrow(pb2022mt)
```


# Complete Results from TWFE Models

## Parameters Estimates Without Pre-matching Processing.

\setcounter{table}{0}
\setcounter{figure}{0}
\def\figurename{Figure E}
\def\tablename{Table E}

Table E\ref{table:2021} and E\ref{table:2022} show the results of the Two-Ways Fixed Effects models for the 2021 presidential election and the 2022 constitutional plebiscite, respectively. Models 1 and 2 of each table contain the unweighted estimates, while the models 3 and 4 show the coefficients from the weighted data. All models are calculated using heteroskedasticity-consistent (robust) standard errors. 

```{r models 2021, results='asis'}

#TWFE OF PRESIDENTIAL ELECTION WITH COMPLETE DATA (WEIGHTED AND NON-WEIGHTED)
m0 <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,data=mm)
m1 <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,data=mm, weights = weight.educ.edad.pre_match)

m0a <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, data=mm)
m1a <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, 
          data=mm, weights = weight.educ.edad.pre_match)

se_m0 <- sqrt(diag(vcov(m0)))
se_m1 <- sqrt(diag(vcovHC(m1,type = "HC0")))
se_m0a <- sqrt(diag(vcov(m0a)))
se_m1a <- sqrt(diag(vcovHC(m1a,type = "HC0")))

pval_m0 <- (1-pnorm(abs(m0$coef/se_m0)))*2
pval_m1 <- (1-pnorm(abs(m1$coef/se_m1)))*2
pval_m0a <- (1-pnorm(abs(m0a$coef/se_m0a)))*2
pval_m1a <- (1-pnorm(abs(m1a$coef/se_m1a)))*2

texreg(l=list(m0,m0a,m1,m1a),
       override.se = list(se_m0,se_m0a,
                          se_m1,se_m1a),
       override.pvalues = list(pval_m0,pval_m0a,
                               pval_m1,pval_m1a),
       custom.coef.names = c("Intercept","Delta Electoral Winner",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models" = 3:4),
       caption = "Pre-Match Weighted and Unweighted Models for 2021 Presidential Election",
       single.row = FALSE, float.pos = "h",use.packages = F,
       label = "table:2021")

```


```{r main effects models, results = 'asis'}
#TWFE OF PLEBISCITE WITH COMPLETE DATA (WEIGHTED AND NON-WEIGHTED)

m0.2 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner, 
                  data=pb2022)
m1.2 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner, 
                  data=pb2022, weights = weight.educ.edad.pre_match)

m0a.2 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01, 
                   data=pb2022)
m1a.2 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01, 
                   data=pb2022, weights = weight.educ.edad.pre_match)

se_m0.2 <- sqrt(diag(vcov(m0.2)))
se_m1.2 <- sqrt(diag(vcovHC(m1.2,type = "HC0")))
se_m0a.2 <- sqrt(diag(vcov(m0a.2)))
se_m1a.2 <- sqrt(diag(vcovHC(m1a.2,type = "HC0")))

pval_m0.2 <- (1-pnorm(abs(m0.2$coef/se_m0.2)))*2
pval_m1.2 <- (1-pnorm(abs(m1.2$coef/se_m1.2)))*2
pval_m0a.2 <- (1-pnorm(abs(m0a.2$coef/se_m0a.2)))*2
pval_m1a.2 <- (1-pnorm(abs(m1a.2$coef/se_m1a.2)))*2

texreg(l=list(m0.2,m0a.2,m1.2,m1a.2), caption = "Pre-Match Weighted and Unweighted Models for 2022 Constitutional Referendum",
       override.se=list(se_m0.2,se_m0a.2,
                        se_m1.2,se_m1a.2),
       override.pvalues = list(pval_m0.2,pval_m0a.2,
                        pval_m1.2,pval_m1a.2),
       custom.coef.names = c("Intercept","Delta Electoral Winner",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       single.row = F, float.pos = "h",use.packages = F,
       label = "table:2022")


```

\newpage

## Parameters Estimates With Pre-matching Processing.

Tables E\ref{table:2021match} and E\ref{table:2022match} show the results from the Two-way Fixed Effects models applied to the from the matched sample of the 2021 presidential and constitutional plebiscite, respectively. In each table Models 1 and 2 contains the unweighted estimates, while the models 3 and 4 show the coefficients from the weighted data. All models are calculated using heteroskedasticity-consistent (robust) standard errors.

```{r weight matching 2021, results='asis'}
#TWFE OF PRESIDENTIAL ELECTION WITH MATCHED DATA (WEIGHTED AND NON-WEIGHTED)

mayo.mt$weights <- NULL

source(paste0(ruta,"3. Analysis Code/Weight_Presidenciales_2021_Match_v2.R"))

mayo.mt <- mayo.mt.w; mayo.mt.w <- NULL

fit0 <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
           data = mayo.mt)
fit1 <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04,
           data = mayo.mt)

fit0a <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
           data = mayo.mt, weights = weight.match.educ.edad)
fit1a <- lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04,
           data = mayo.mt, weights = weight.match.educ.edad)

se_fit0 <- sqrt(diag(vcov(fit0)))
se_fit1 <- sqrt(diag(vcov(fit1)))

se_fit0a <- sqrt(diag(vcovHC(fit0a,type="HC0")))
se_fit1a <- sqrt(diag(vcovHC(fit1a,type="HC0")))

pval_fit0 <- (1-pnorm(abs(fit0$coef/se_fit0)))*2
pval_fit1 <- (1-pnorm(abs(fit1$coef/se_fit1)))*2

pval_fit0a <- (1-pnorm(abs(fit0a$coef/se_fit0a)))*2
pval_fit1a <- (1-pnorm(abs(fit1a$coef/se_fit1a)))*2

texreg(list(fit0,fit1,fit0a,fit1a),digitis=4,custom.coef.names = c("Intercept","Delta Electoral Winner",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       override.se = list(se_fit0,se_fit1,se_fit0a,se_fit1a),
       override.pvalues = list(pval_fit0,pval_fit1,pval_fit0a,pval_fit1a),
       custom.header=list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       caption = "Matched Weighted and Unweighted Models for 2021 Presidential Election",
       single.row = F, float.pos = "h",label = "table:2021match",use.packages = F)
```

```{r matching model,results='asis'}
#TWFE OF PLEBISCITE WITH MATCHED DATA (WEIGHTED AND NON-WEIGHTED)

pb2022mt$weights <- NULL

source(paste0(ruta,"3. Analysis Code/Weight_Plebiscito_Match_v2.R"))

pb2022mt <- pb2022mt.w; pb2022mt.w <- NULL

#Models
fit3 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,
           data = pb2022mt)
fit4 <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01,
           data = pb2022mt)

fit3a <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,
           data = pb2022mt, weights = weight.match.educ.edad)
fit4a <- lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01,
           data = pb2022mt, weights = weight.match.educ.edad)

se_fit3 <- sqrt(diag(vcov(fit3)))
se_fit4 <- sqrt(diag(vcov(fit4)))

se_fit3a <- sqrt(diag(vcovHC(fit3a,type="HC0")))
se_fit4a <- sqrt(diag(vcovHC(fit4a,type="HC0")))

pval_fit3 <- (1-pnorm(abs(fit3$coef/se_fit3)))*2
pval_fit4 <- (1-pnorm(abs(fit4$coef/se_fit4)))*2

pval_fit3a <- (1-pnorm(abs(fit3a$coef/se_fit3a)))*2
pval_fit4a <- (1-pnorm(abs(fit4a$coef/se_fit4a)))*2



texreg(l=list(fit3,fit4,fit3a,fit4a), caption = "Matched Weighted and Unweighted Models for 2022 Constitutional Referendum",
       override.se = list(se_fit3,se_fit4,
                          se_fit3a,se_fit4a),
       override.pvalues = list(pval_fit3,pval_fit4,
                               pval_fit3a,pval_fit4a),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       custom.coef.names =  c("Intercept","Delta Electoral Winner",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       single.row = F, float.pos = "h",use.packages = F,label = "table:2022match")

```

\newpage


# Auto-regressive models

\setcounter{table}{0}
\setcounter{figure}{0}
\def\figurename{Figure F}
\def\tablename{Table F}


## Parameters Estimates Without Pre-matching Processing.


```{r auto-regressive models,results='asis'}

#AUTOREGRESSIVE MODEL OF PRESIDENTIAL ELECTION WITH COMPLETE DATA (WEIGHTED AND NON-WEIGHTED)
m0 <- lm(conf_soc_w05 ~ delta_elec_winner+conf_soc_w04,data=mm)
m1 <- lm(conf_soc_w05 ~ delta_elec_winner+conf_soc_w04,data=mm, weights = weight.educ.edad.pre_match)

m0a <- lm(conf_soc_w05 ~ delta_elec_winner:pol_interesF2_w04+conf_soc_w04, data=mm)
m1a <- lm(conf_soc_w05 ~ delta_elec_winner:pol_interesF2_w04+conf_soc_w04, 
          data=mm, weights = weight.educ.edad.pre_match)

se_m0 <- sqrt(diag(vcov(m0)))
se_m1 <- sqrt(diag(vcovHC(m1,type = "HC0")))
se_m0a <- sqrt(diag(vcov(m0a)))
se_m1a <- sqrt(diag(vcovHC(m1a,type = "HC0")))

pval_m0 <- (1-pnorm(abs(m0$coef/se_m0)))*2
pval_m1 <- (1-pnorm(abs(m1$coef/se_m1)))*2
pval_m0a <- (1-pnorm(abs(m0a$coef/se_m0a)))*2
pval_m1a <- (1-pnorm(abs(m1a$coef/se_m1a)))*2

texreg(l=list(m0,m0a,m1,m1a),
       override.se = list(se_m0,se_m0a,
                          se_m1,se_m1a),
       override.pvalues = list(pval_m0,pval_m0a,
                               pval_m1,pval_m1a),
       custom.coef.names = c("Intercept","Delta Electoral Winner","Social Trust (t-1)",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models" = 3:4),
       caption = "Pre-Match Weighted and Unweighted Auto-regressive Models for 2021 Presidential Election",
       single.row = FALSE, float.pos = "h!",use.packages = F,
       label = "table:auto2021")

```


```{r,results='asis'}
#AUTOREGRESSIVE MODEL OF PLEBISCITE WITH COMPLETE DATA (WEIGHTED AND NON-WEIGHTED)

m0.2 <- lm(prs_trst_w02 ~ delta_elec_winner + prs_trst_w01, 
                  data=pb2022)
m1.2 <- lm(prs_trst_w02 ~ delta_elec_winner + prs_trst_w01, 
                  data=pb2022, weights = weight.educ.edad.pre_match)

m0a.2 <- lm(prs_trst_w02 ~ delta_elec_winner:pol_interesF2_w01 + prs_trst_w01, 
                   data=pb2022)
m1a.2 <- lm(prs_trst_w02 ~ delta_elec_winner:pol_interesF2_w01 + prs_trst_w01, 
                   data=pb2022, weights = weight.educ.edad.pre_match)

se_m0.2 <- sqrt(diag(vcov(m0.2)))
se_m1.2 <- sqrt(diag(vcovHC(m1.2,type = "HC0")))
se_m0a.2 <- sqrt(diag(vcov(m0a.2)))
se_m1a.2 <- sqrt(diag(vcovHC(m1a.2,type = "HC0")))

pval_m0.2 <- (1-pnorm(abs(m0.2$coef/se_m0.2)))*2
pval_m1.2 <- (1-pnorm(abs(m1.2$coef/se_m1.2)))*2
pval_m0a.2 <- (1-pnorm(abs(m0a.2$coef/se_m0a.2)))*2
pval_m1a.2 <- (1-pnorm(abs(m1a.2$coef/se_m1a.2)))*2

texreg(l=list(m0.2,m0a.2,m1.2,m1a.2), caption = "Pre-Match Weighted and Unweighted Models for 2022 Constitutional Referendum",
       override.se=list(se_m0.2,se_m0a.2,
                        se_m1.2,se_m1a.2),
       override.pvalues = list(pval_m0.2,pval_m0a.2,
                        pval_m1.2,pval_m1a.2),
       custom.coef.names = c("Intercept","Delta Electoral Winner","Social Trust (t-1)",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       single.row = F, float.pos = "h!",use.packages = F,
       label = "table:auto2022")

```
\newpage

## Parameters Estimates With Pre-matching Processing.

```{r,results='asis'}
#AUTOREGRESSIVE MODEL OF PRESIDENTIAL ELECTION WITH MATCHED DATA (WEIGHTED AND NON-WEIGHTED)

fit0 <- lm(conf_soc_w05 ~ delta_elec_winner+conf_soc_w04,
           data = mayo.mt)
fit1 <- lm(conf_soc_w05 ~ delta_elec_winner:pol_interesF2_w04+conf_soc_w04,
           data = mayo.mt)

fit0a <- lm(conf_soc_w05 ~ delta_elec_winner+conf_soc_w04,
           data = mayo.mt, weights = weight.match.educ.edad)
fit1a <- lm(conf_soc_w05 ~ delta_elec_winner:pol_interesF2_w04+conf_soc_w04,
           data = mayo.mt, weights = weight.match.educ.edad)

se_fit0 <- sqrt(diag(vcov(fit0)))
se_fit1 <- sqrt(diag(vcov(fit1)))

se_fit0a <- sqrt(diag(vcovHC(fit0a,type="HC0")))
se_fit1a <- sqrt(diag(vcovHC(fit1a,type="HC0")))

pval_fit0 <- (1-pnorm(abs(fit0$coef/se_fit0)))*2
pval_fit1 <- (1-pnorm(abs(fit1$coef/se_fit1)))*2

pval_fit0a <- (1-pnorm(abs(fit0a$coef/se_fit0a)))*2
pval_fit1a <- (1-pnorm(abs(fit1a$coef/se_fit1a)))*2

texreg(list(fit0,fit1,fit0a,fit1a),digitis=4,
       custom.coef.names =c("Intercept","Delta Electoral Winner","Social Trust (t-1)",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       override.se = list(se_fit0,se_fit1,se_fit0a,se_fit1a),
       override.pvalues = list(pval_fit0,pval_fit1,pval_fit0a,pval_fit1a),
       custom.header=list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       caption = "Matched Weighted and Unweighted Auto-regressive Models for 2021 Presidential Election",
       single.row = F, float.pos = "h!",label = "table:auto2021match",use.packages = F)
```



```{r, results='asis'}

#AUTOREGRESSIVE MODEL OF PLEBISCITE WITH MATCHED DATA (WEIGHTED AND NON-WEIGHTED)

fit3 <- lm(prs_trst_w02 ~ delta_elec_winner+prs_trst_w01,
           data = pb2022mt)
fit4 <- lm(prs_trst_w02 ~ delta_elec_winner:pol_interesF2_w01+prs_trst_w01,
           data = pb2022mt)

fit3a <- lm(prs_trst_w02 ~ delta_elec_winner+prs_trst_w01,
           data = pb2022mt, weights = weight.match.educ.edad)
fit4a <- lm(prs_trst_w02 ~ delta_elec_winner:pol_interesF2_w01+prs_trst_w01,
           data = pb2022mt, weights = weight.match.educ.edad)

se_fit3 <- sqrt(diag(vcov(fit3)))
se_fit4 <- sqrt(diag(vcov(fit4)))

se_fit3a <- sqrt(diag(vcovHC(fit3a,type="HC0")))
se_fit4a <- sqrt(diag(vcovHC(fit4a,type="HC0")))

pval_fit3 <- (1-pnorm(abs(fit3$coef/se_fit3)))*2
pval_fit4 <- (1-pnorm(abs(fit4$coef/se_fit4)))*2

pval_fit3a <- (1-pnorm(abs(fit3a$coef/se_fit3a)))*2
pval_fit4a <- (1-pnorm(abs(fit4a$coef/se_fit4a)))*2



texreg(l=list(fit3,fit4,fit3a,fit4a), 
       caption = "Matched Weighted and Unweighted Auto-regressive Models for 2022 Constitutional Referendum",
       override.se = list(se_fit3,se_fit4,
                          se_fit3a,se_fit4a),
       override.pvalues = list(pval_fit3,pval_fit4,
                               pval_fit3a,pval_fit4a),
       custom.header = list("Unweighted Models"=1:2,"Weighted Models"=3:4),
       custom.coef.names = c("Intercept","Delta Electoral Winner","Social Trust (t-1)",
                             "Delta Electoral Winner*Low Political Interest",
                             "Delta Electoral Winner*Mid Political Interest",
                             "Delta Electoral Winner*High Political Interest"),
       single.row = F, float.pos = "h!",use.packages = F,label = "table:auto2022match")
```





\newpage

# TWFE parameters with different tolerance levels

\setcounter{table}{0}
\setcounter{figure}{0}
\def\figurename{Figure G}
\def\tablename{Table G}


```{r}

### Loop para match: se cambian los umbrales de tolerancia, 0.01 a 0.5

#####---- PRESIDENTIAL -----####

results_loops <- list()
regressions <- list()
match_df <- list()
love_plots <- list()
boric <- list()
pleb <- list()

fits_boric <- list()

tolerance <- c(0.01,0.05,0.10,0.20,0.30,0.40,0.50)

# Eleccion de Boric

for (i in 1:length(tolerance)){
  tols <- tolerance[i]
  
  m1.full <- matchit(delta_elec_winner ~ lrscaleF_w04 + educaF_w03 + edad_w04F + mujer_w03 + 
                       pol_interesF_w04, data=mm, method = "cardinality", tols = tols, solver="gurobi")
  
  v <- data.frame(old = c("lrscaleF_w04_Izquierda","lrscaleF_w04_Centro","lrscaleF_w04_Derecha",
                          "lrscaleF_w04_No Identificado","educaF_w03_Secondary Ed.",
                          "educaF_w03_Tertiary Tech.",
                          "educaF_w03_Tertiary Univ","edad_w04F_18-29","edad_w04F_30-39",
                          "edad_w04F_40-49","edad_w04F_50-65",
                          "mujer_w03","pol_interesF_w04","distance"),
                  new = c("Left", "Center", "Right", "No Political Id.", 
                          "Secondary Ed.", "Tertiary Tech.", "Tertiary Univ.", "18-29","30-39","40-49","50-65",
                          "Women","No Political Interest","Propensity Score"))
  
  love_plot <- love.plot(m1.full, abs=T, title="Covariate Balance for Electoral Winning",
            sample.names=c("Unmatched", "Matched"), stars = "std",
            var.names = v)
  
  mayo.mt <- match.data(m1.full)
  row.names(mayo.mt) <- 1:nrow(mayo.mt)
  mayo.mt$weights <- NULL
  
  source(paste0(ruta, "3. Analysis Code/Weight_Presidenciales_2021_Match_v2.R"))
  
  mayo.mt <- mayo.mt.w; mayo.mt.w <- NULL
  
  #No ponderado
  
  fit1 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
                        data = mayo.mt),vcov=vcovHC)
  fit2 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF_w04,
                        data = mayo.mt),vcov=vcovHC)
  
  regr <- screenreg(l=list(fit1,fit2),digits=3,single_row=F,
                    custom.model.names = c("Non-Weighted Model 1","Non-Weighted Model 2"),
                    stars=c(0.001,0.01,0.05,0.1))
  
  ## Ponderado
  
  fit1_w <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
             data = mayo.mt, weights = weight.match.educ.edad),vcov=vcovHC)
  fit2_w <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF_w04,
             data = mayo.mt, weights = weight.match.educ.edad),vcov=vcovHC)
  
  regr_w <- screenreg(l=list(fit1_w,fit2_w),digits=3,single_row=F,
                      custom.model.names = c("Weighted Model 1","Weighted Model 2"),
                      stars=c(0.001,0.01,0.05,0.1))
  
  # Storeage
  
  love_plots[[i]] <- love_plot
  names(love_plots)[i] <- paste("tolerance:",tols,sep=" ")
  
  match_df[[i]] <- mayo.mt
  names(match_df)[i] <- paste("tolerance:",tols,sep=" ")
  
  regressions[[i]] <- setNames(list(regr,regr_w),c("Non weighted model","Weighted model"))
  names(regressions)[i] <- paste("tolerance:",tols,sep=" ")
  
  boric[[i]] <- setNames(list(love_plots,match_df,regressions),
                                 c("Love Plots","Matched DF","Regression"))
  names(boric)[i] <- paste("tolerance:",tols,sep=" ")
  
  ## Para el screenreg de todas las regresiones
  
  fits_boric[[i]] <- list(fit1,fit2,fit1_w,fit2_w)
}


t.01 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.01"]]$lrscaleF_w04))*100)
t.01$tolerance = 0.01

t.05 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.05"]]$lrscaleF_w04))*100)
t.05$tolerance = 0.05

t.1 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.1"]]$lrscaleF_w04))*100)
t.1$tolerance <- 0.1

t.2 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.2"]]$lrscaleF_w04))*100)
t.2$tolerance <- 0.2

t.3 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.3"]]$lrscaleF_w04))*100)
t.3$tolerance <- 0.3

t.4 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.4"]]$lrscaleF_w04))*100)
t.4$tolerance <- 0.4

t.5 <- data.frame(prop.table(table(boric[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.5"]]$lrscaleF_w04))*100)
t.5$tolerance <- 0.5

tolerance.boric <- rbind(t.01,t.05,t.1,t.2,t.3,t.4,t.5)


#ggplot(tolerance.boric,aes(x=Var1,y=Freq)) + geom_bar(stat="identity") + facet_wrap(~tolerance) + 
#  labs(title = "Presidential Election",x="Political Identification") +
#  scale_x_discrete(labels=c("Izquierda"="Left",
#                          "Centro"="Center",
#                          "Derecha"="Right",
#                          "No Identificado"="Independent"))


## Unweighted Estimates

#screenreg(list(fits_boric[[1]][[1]],fits_boric[[1]][[2]],
#               fits_boric[[2]][[1]],fits_boric[[2]][[2]],
#               fits_boric[[3]][[1]],fits_boric[[3]][[2]],
#               fits_boric[[4]][[1]],fits_boric[[4]][[2]],
#               fits_boric[[5]][[1]],fits_boric[[5]][[2]],
#               fits_boric[[6]][[1]],fits_boric[[6]][[2]],
#               fits_boric[[7]][[1]],fits_boric[[7]][[2]]),
#          digits = 3, custom.header = list("Tolerance = 0.01" = 1:2,
#                                    "Tolerance = 0.05" = 3:4,
#                                    "Tolerance = 0.1" = 5:6,
#                                    "Tolerance = 0.2" = 7:8,
#                                    "Tolerance = 0.3" = 9:10,
#                                    "Tolerance = 0.4" = 11:12,
#                                    "Tolerance = 0.5" = 13:14))

# Estimaciones no ponderadas

coef.df.list <- list()

for (i in 1:7) {
  coef.df <- data.frame(coef = fits_boric[[i]][[1]][,1],
                        se = fits_boric[[i]][[1]][,2],
                        pr = fits_boric[[i]][[1]][,4])
  
  coef.df$stars <- ifelse(coef.df$pr < 0.001, "***",
                          ifelse(coef.df$pr < 0.01, "**",
                                 ifelse(coef.df$pr < 0.05, "*", "")))
  
  coef.df$n = nobs(fits_boric[[i]][[1]])
  
  coef.df$tolerance = rep(tolerance[i])
  
  t.boric <- tolerance.boric %>% filter(tolerance==coef.df$tolerance)
  
  coef.df$left <- NA
  coef.df$center <- NA
  coef.df$right <- NA
  coef.df$non.ideologues <- NA
  
  coef.df$left[2] <- t.boric[1,2]
  coef.df$center[2] <- t.boric[2,2]
  coef.df$right[2] <- t.boric[3,2]
  coef.df$non.ideologues[2] <- t.boric[4,2]
  
  coef.df.list[[i]] <- coef.df
}

# Combine the data frames
combined_df <- do.call(rbind, coef.df.list)
combined_df <- combined_df %>% rownames_to_column(var="Estimate")

# Resultados sin Match

m1.sm <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
                     data = mm),vcov=vcovHC)

coef.df <- data.frame(coef = m1.sm[,1],
                      se = m1.sm[,2],
                      pr = m1.sm[,4])

coef.df$stars <- ifelse(coef.df$pr < 0.001, "***",
                        ifelse(coef.df$pr < 0.01, "**",
                               ifelse(coef.df$pr < 0.05, "*", "")))

coef.df$n = nobs(m1.sm)

coef.df$tolerance = "No Match"

coef.df$left <- NA
coef.df$center <- NA
coef.df$right <- NA
coef.df$non.ideologues <- NA

coef.df$left[2] <- data.frame(prop.table(table(mm$lrscaleF_w04))*100)[1,2]
coef.df$center[2] <- data.frame(prop.table(table(mm$lrscaleF_w04))*100)[2,2]
coef.df$right[2] <- data.frame(prop.table(table(mm$lrscaleF_w04))*100)[3,2]
coef.df$non.ideologues[2] <- data.frame(prop.table(table(mm$lrscaleF_w04))*100)[4,2]

coef.df <- coef.df %>% rownames_to_column(var="Estimate") %>% 
  mutate(Estimate = paste(Estimate,"7",sep=""))

combined_df <- rbind(combined_df,coef.df)
combined_df <- combined_df %>% filter(!Estimate %in% c("(Intercept)","(Intercept)1","(Intercept)2","(Intercept)3",
                                                    "(Intercept)4","(Intercept)5","(Intercept)6","(Intercept)7")) %>%
  mutate(Estimate = "Electoral status") %>% select(!stars)

cols <- c("","Coef","SE","p-value","n","tol","Left %","Center %","Right %","Non-ideo %")

colnames(combined_df) <- cols

print(xtable(combined_df, digits=c(0, 0, rep(2, length(colnames(combined_df))-1)),
             align=c("l","c","c","c", "p{1.5cm}","p{1.5cm}","c","c","c","c","c"),
             caption = "Point estimates of winner-loser gap on social trust at different tolerance levels for the Presidential Election with Unweighted data",
             label = "table:tolerance_estimate"),include.rownames=FALSE,
      scalebox=.85, table.placement="H",comment = FALSE)



#####---- Plebiscito -----####

regressions <- list()
match_df <- list()
love_plots <- list()

fits_plbs <- list()

for (i in 1:length(tolerance)){
  tols <- tolerance[i]
  
  m1.full <- matchit(delta_elec_winner ~ lrscaleF_w01 +
                       educaF_w01 + edad_w01F + mujer_w01 + pol_interesF_w01, 
                     data=subset(pb2022, !is.na(delta_elec_winner)), 
                     #method = "full", replace=T)
                     method = "cardinality", tols = tols, solver="gurobi")
  
  v <- data.frame(old = c("lrscaleF_w01_Izquierda","lrscaleF_w01_Centro","lrscaleF_w01_Derecha",
                          "lrscaleF_w01_No Identificado","educaF_w01_Secondary Ed.",
                          "educaF_w01_Tertiary Tech.",
                          "educaF_w01_Tertiary Univ","edad_w01F_18-29","edad_w01F_30-39",
                          "edad_w01F_40-49","edad_w01F_50-65",
                          "mujer_w01","pol_interesF_w01","distance"),
                  new = c("Left", "Center", "Right", "No Political Id.", 
                          "Secondary Ed.", "Tertiary Tech.", "Tertiary Univ.", "18-29","30-39","40-49","50-65",
                          "Women","No Political Interest","Propensity Score"))
  
  love_plot <- love.plot(m1.full, abs=T, 
                         title="Covariate Balance for Electoral Winning",
                         sample.names=c("Unmatched","Matched"))
  
  pb2022mt <- match.data(m1.full)
  pb2022mt$weights <- NULL
  
  source(paste0(ruta, "3. Analysis Code/Weight_Plebiscito_Match_v2.R"))

  pb2022mt <- pb2022mt.w; pb2022mt.w <- NULL
  
  #No ponderado
  
  fit1 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,
                        data = pb2022mt),vcov=vcovHC)
  fit2 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF_w01,
                        data = pb2022mt),vcov=vcovHC)
  
  regr <- screenreg(l=list(fit1,fit2),digits=3,single_row=F,
                    custom.model.names = c("Non-Weighted Model 1","Non-Weighted Model 2"),
                    stars=c(0.001,0.01,0.05,0.1))
  
  ## Ponderado
  
  fit1_w <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,
                        data = pb2022mt, weights = weight.match.educ.edad),vcov=vcovHC)
  fit2_w <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF_w01,
                        data = pb2022mt, weights = weight.match.educ.edad),vcov=vcovHC)
  
  regr_w <- screenreg(l=list(fit1_w,fit2_w),digits=3,single_row=F,
                      custom.model.names = c("Weighted Model 1","Weighted Model 2"),
                      stars=c(0.001,0.01,0.05,0.1))
  
  # Storeage
  
  love_plots[[i]] <- love_plot
  names(love_plots)[i] <- paste("tolerance:",tols,sep=" ")
  
  match_df[[i]] <- pb2022mt
  names(match_df)[i] <- paste("tolerance:",tols,sep=" ")
  
  regressions[[i]] <- setNames(list(regr,regr_w),c("Non weighted model","Weighted model"))
  names(regressions)[i] <- paste("tolerance:",tols,sep=" ")
  
  pleb[[i]] <- setNames(list(love_plots,match_df,regressions),
                                 c("Love Plots","Matched DF","Regression"))
  names(pleb)[i] <- paste("tolerance:",tols,sep=" ")
  
  fits_plbs[[i]] <- list(fit1,fit2,fit1_w,fit2_w)
}

## Unweighted Estimates

t.01 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.01"]]$lrscaleF_w01))*100)
t.01$tolerance = 0.01

t.05 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.05"]]$lrscaleF_w01))*100)
t.05$tolerance = 0.05

t.1 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.1"]]$lrscaleF_w01))*100)
t.1$tolerance <- 0.1

t.2 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.2"]]$lrscaleF_w01))*100)
t.2$tolerance <- 0.2

t.3 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.3"]]$lrscaleF_w01))*100)
t.3$tolerance <- 0.3

t.4 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.4"]]$lrscaleF_w01))*100)
t.4$tolerance <- 0.4

t.5 <- data.frame(prop.table(table(pleb[["tolerance: 0.5"]][["Matched DF"]][["tolerance: 0.5"]]$lrscaleF_w01))*100)
t.5$tolerance <- 0.5

tolerance.pleb <- rbind(t.01,t.05,t.1,t.2,t.3,t.4,t.5)
coef.df.list <- list()

for (i in 1:7) {
  coef.df <- data.frame(coef = fits_plbs[[i]][[1]][,1],
                        se = fits_plbs[[i]][[1]][,2],
                        pr = fits_plbs[[i]][[1]][,4])
  
  coef.df$stars <- ifelse(coef.df$pr < 0.001, "***",
                          ifelse(coef.df$pr < 0.01, "**",
                                 ifelse(coef.df$pr < 0.05, "*", "")))
  
  coef.df$n = nobs(fits_plbs[[i]][[1]])
  
  coef.df$tolerance = rep(tolerance[i])
  
  t.pleb <- tolerance.pleb %>% filter(tolerance==coef.df$tolerance)
  
  coef.df$left <- NA
  coef.df$center <- NA
  coef.df$right <- NA
  coef.df$non.ideologues <- NA
  
  coef.df$left[2] <- t.pleb[1,2]
  coef.df$center[2] <- t.pleb[2,2]
  coef.df$right[2] <- t.pleb[3,2]
  coef.df$non.ideologues[2] <- t.pleb[4,2]
  
  coef.df.list[[i]] <- coef.df
}

# Combine the data frames
combined_df <- do.call(rbind, coef.df.list)
combined_df <- combined_df %>% rownames_to_column(var="Estimate")

# Resultados sin Match

m1.sm <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,
                     data = pb2022),vcov=vcovHC)

coef.df <- data.frame(coef = m1.sm[,1],
                      se = m1.sm[,2],
                      pr = m1.sm[,4])

coef.df$stars <- ifelse(coef.df$pr < 0.001, "***",
                        ifelse(coef.df$pr < 0.01, "**",
                               ifelse(coef.df$pr < 0.05, "*", "")))

coef.df$n = nobs(m1.sm)

coef.df$tolerance = "No Match"

coef.df$left <- NA
coef.df$center <- NA
coef.df$right <- NA
coef.df$non.ideologues <- NA

coef.df$left[2] <- data.frame(prop.table(table(pb2022$lrscaleF_w01))*100)[1,2]
coef.df$center[2] <- data.frame(prop.table(table(pb2022$lrscaleF_w01))*100)[2,2]
coef.df$right[2] <- data.frame(prop.table(table(pb2022$lrscaleF_w01))*100)[3,2]
coef.df$non.ideologues[2] <- data.frame(prop.table(table(pb2022$lrscaleF_w01))*100)[4,2]

coef.df <- coef.df %>% rownames_to_column(var="Estimate") %>% 
  mutate(Estimate = paste(Estimate,"7",sep=""))

combined_df <- rbind(combined_df,coef.df)
combined_df <- combined_df %>% filter(!Estimate %in% c("(Intercept)","(Intercept)1","(Intercept)2","(Intercept)3",
"(Intercept)4","(Intercept)5","(Intercept)6","(Intercept)7")) %>%
  mutate(Estimate = "Electoral status") %>% select(!stars)

cols <- c("","Coef","SE","p-value","n","tol","Left %","Center %","Right %","Non-ideo %")

colnames(combined_df) <- cols

print(xtable(combined_df, digits=c(0, 0, rep(2, length(colnames(combined_df))-1)),
             align=c("l","c","c","c", "p{1.5cm}","p{1.5cm}","c","c","c","c","c"),
             caption = "Point estimates of winner-loser gap on social trust at different tolerance levels for the Constitutional Plebiscite with Unweighted data",
             label = "table:tolerance_estimate"),include.rownames = FALSE,
      scalebox=.85, table.placement="H",comment = FALSE)

#save(results_loops, file=paste0(ruta,"listas.RData"))

```

\newpage

# Replication of Figure 2


```{r graficos modelos, eval=FALSE, fig.cap="\\label{fig:effects}Replication of Figure 2", fig.height=7, fig.pos="h", fig.width=10, include=FALSE, results='asis'}

### Pres

## Main effect

fit0 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner, data=mm,
                    weights = weight.educ.edad.pre_match),
                 vcov=vcovHC)
vars <- c("Intercept","Electoral Winner")
coef <- fit0[,1]
se <- fit0[,2]
df <- data.frame(cbind(vars,coef,se))
df$cimin95 <- NA
df$cimin95 <- as.numeric(df$coef)-1.96*as.numeric(df$se)
df$cimax95 <- NA
df$cimax95 <- as.numeric(df$coef)+1.96*as.numeric(df$se)
df$cimin90 <- NA
df$cimin90 <- as.numeric(df$coef)-1.65*as.numeric(df$se)
df$cimax90 <- NA
df$cimax90 <- as.numeric(df$coef)+1.65*as.numeric(df$se)
df$model <- "Weighted Unmatched \nSample"
df$n <- nobs(fit0)

fit1 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
         data = mayo.mt, weights = weight.match.educ.edad),vcov=vcovHC)

coef_2 <- fit1[,1]
se_2 <- fit1[,2]
df_2 <- data.frame(cbind(vars,coef=coef_2,se=se_2))
df_2$cimin95 <- NA
df_2$cimin95 <- as.numeric(df_2$coef)-1.96*as.numeric(df_2$se)
df_2$cimax95 <- NA
df_2$cimax95 <- as.numeric(df_2$coef)+1.96*as.numeric(df_2$se)
df_2$cimin90 <- NA
df_2$cimin90 <- as.numeric(df_2$coef)-1.65*as.numeric(df_2$se)
df_2$cimax90 <- NA
df_2$cimax90 <- as.numeric(df_2$coef)+1.65*as.numeric(df_2$se)
df_2$model <- "Weighted Matched \nSample"
df_2$n <- nobs(fit1)

#############################################################################
#############################################################################


fit0_uw <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner, data=mm),
                    vcov=vcovHC)

vars <- c("Intercept","Electoral Winner")
coef <- fit0_uw[,1]
se <- fit0_uw[,2]
df_uw <- data.frame(cbind(vars,coef,se))
df_uw$cimin95 <- NA
df_uw$cimin95 <- as.numeric(df_uw$coef)-1.96*as.numeric(df_uw$se)
df_uw$cimax95 <- NA
df_uw$cimax95 <- as.numeric(df_uw$coef)+1.96*as.numeric(df_uw$se)
df_uw$cimin90 <- NA
df_uw$cimin90 <- as.numeric(df_uw$coef)-1.65*as.numeric(df_uw$se)
df_uw$cimax90 <- NA
df_uw$cimax90 <- as.numeric(df_uw$coef)+1.65*as.numeric(df_uw$se)
df_uw$model <- "Unweighted Unmatched \nSample"
df_uw$n <- nobs(fit0_uw)


fit1_uw <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner,
                    data = mayo.mt),vcov=vcovHC)


coef_2 <- fit1_uw[,1]
se_2 <- fit1_uw[,2]
df_uw_2 <- data.frame(cbind(vars,coef=coef_2,se=se_2))
df_uw_2$cimin95 <- NA
df_uw_2$cimin95 <- as.numeric(df_uw_2$coef)-1.96*as.numeric(df_uw_2$se)
df_uw_2$cimax95 <- NA
df_uw_2$cimax95 <- as.numeric(df_uw_2$coef)+1.96*as.numeric(df_uw_2$se)
df_uw_2$cimin90 <- NA
df_uw_2$cimin90 <- as.numeric(df_uw_2$coef)-1.65*as.numeric(df_uw_2$se)
df_uw_2$cimax90 <- NA
df_uw_2$cimax90 <- as.numeric(df_uw_2$coef)+1.65*as.numeric(df_uw_2$se)
df_uw_2$model <- "Unweighted Matched \nSample"
df_uw_2$n <- nobs(fit1_uw)

df_a_w <- rbind(df,df_2)
df_a_uw <- rbind(df_uw,df_uw_2)

## Interaction with political interest

fit2 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, data=mm,
                    weights = weight.educ.edad.pre_match),vcov=vcovHC)



vars <- c("Intercept","Low Interest","Medium Interest","High Interest")
coef <- fit2[,1]
se <- fit2[,2]
df <- data.frame(cbind(vars,coef,se))
df$cimin95 <- NA
df$cimin95 <- as.numeric(df$coef)-1.96*as.numeric(df$se)
df$cimax95 <- NA
df$cimax95 <- as.numeric(df$coef)+1.96*as.numeric(df$se)
df$cimin90 <- NA
df$cimin90 <- as.numeric(df$coef)-1.65*as.numeric(df$se)
df$cimax90 <- NA
df$cimax90 <- as.numeric(df$coef)+1.65*as.numeric(df$se)
df$model <- "Weighted Unmatched \nSample"
df$n <- nobs(fit2)

fit3 <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, data = mayo.mt, 
                    weights = weight.match.educ.edad),
                 vcov=vcovHC)


coef <- fit3[,1]
se <- fit3[,2]
df_2 <- data.frame(cbind(vars,coef,se))
df_2$cimin95 <- NA
df_2$cimin95 <- as.numeric(df_2$coef)-1.96*as.numeric(df_2$se)
df_2$cimax95 <- NA
df_2$cimax95 <- as.numeric(df_2$coef)+1.96*as.numeric(df_2$se)
df_2$cimin90 <- NA
df_2$cimin90 <- as.numeric(df_2$coef)-1.65*as.numeric(df_2$se)
df_2$cimax90 <- NA
df_2$cimax90 <- as.numeric(df_2$coef)+1.65*as.numeric(df_2$se)
df_2$model <- "Weighted Matched \nSample"
df_2$n <- nobs(fit3)


fit2_uw <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, data=mm),
                    vcov=vcovHC)



vars <- c("Intercept","Low Interest","Medium Interest","High Interest")
coef <- fit2_uw[,1]
se <- fit2_uw[,2]
df_uw <- data.frame(cbind(vars,coef,se))
df_uw$cimin95 <- NA
df_uw$cimin95 <- as.numeric(df_uw$coef)-1.96*as.numeric(df_uw$se)
df_uw$cimax95 <- NA
df_uw$cimax95 <- as.numeric(df_uw$coef)+1.96*as.numeric(df_uw$se)
df_uw$cimin90 <- NA
df_uw$cimin90 <- as.numeric(df_uw
                            $coef)-1.65*as.numeric(df_uw$se)
df_uw$cimax90 <- NA
df_uw$cimax90 <- as.numeric(df_uw$coef)+1.65*as.numeric(df_uw$se)
df_uw$model <- "Unweighted Unmatched \nSample"
df_uw$n <- nobs(fit2_uw)

fit3_uw <- coeftest(lm(conf_soc_w05-conf_soc_w04 ~ delta_elec_winner:pol_interesF2_w04, data = mayo.mt),
                 vcov=vcovHC)


coef <- fit3_uw[,1]
se <- fit3_uw[,2]
df_uw_2 <- data.frame(cbind(vars,coef,se))
df_uw_2$cimin95 <- NA
df_uw_2$cimin95 <- as.numeric(df_uw_2$coef)-1.96*as.numeric(df_uw_2$se)
df_uw_2$cimax95 <- NA
df_uw_2$cimax95 <- as.numeric(df_uw_2$coef)+1.96*as.numeric(df_uw_2$se)
df_uw_2$cimin90 <- NA
df_uw_2$cimin90 <- as.numeric(df_uw_2$coef)-1.65*as.numeric(df_uw_2$se)
df_uw_2$cimax90 <- NA
df_uw_2$cimax90 <- as.numeric(df_uw_2$coef)+1.65*as.numeric(df_uw_2$se)
df_uw_2$model <- "Unweighted Matched \nSample"
df_uw_2$n <- nobs(fit3_uw)


df_b_w <- rbind(df,df_2)
df_b_uw <- rbind(df_uw,df_uw_2)

df_b_w$effect <- "Interactive Effects"
df_b_uw$effect <- "Interactive Effects"

df_a_w$effect <- "Main Effects"
df_a_uw$effect <- "Main Effects"

df_pr_uw <- rbind(df_a_uw,df_b_uw) # Unweighted estimates
df_pr_w <- rbind(df_a_w,df_b_w) # Weighted estimates



### Pleb

fit0 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,data=pb2022,
                    weights = weight.educ.edad.pre_match),vocv=vcovHC)

vars <- c("Intercept","Electoral Winner")
coef <- fit0[,1]
se <- fit0[,2]
df <- data.frame(cbind(vars,coef,se))
df$cimin95 <- NA
df$cimin95 <- as.numeric(df$coef)-1.96*as.numeric(df$se)
df$cimax95 <- NA
df$cimax95 <- as.numeric(df$coef)+1.96*as.numeric(df$se)
df$cimin90 <- NA
df$cimin90 <- as.numeric(df$coef)-1.65*as.numeric(df$se)
df$cimax90 <- NA
df$cimax90 <- as.numeric(df$coef)+1.65*as.numeric(df$se)
df$model <- "Weighted Unmatched \nSample"
df$n <- nobs(fit0)


fit1 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,data=pb2022mt,
                    weights = weight.match.educ.edad),vcov=vcovHC)

coef_2 <- fit1[,1]
se_2 <- fit1[,2]
df_2 <- data.frame(cbind(vars,coef=coef_2,se=se_2))
df_2$cimin95 <- NA
df_2$cimin95 <- as.numeric(df_2$coef)-1.96*as.numeric(df_2$se)
df_2$cimax95 <- NA
df_2$cimax95 <- as.numeric(df_2$coef)+1.96*as.numeric(df_2$se)
df_2$cimin90 <- NA
df_2$cimin90 <- as.numeric(df_2$coef)-1.65*as.numeric(df_2$se)
df_2$cimax90 <- NA
df_2$cimax90 <- as.numeric(df_2$coef)+1.65*as.numeric(df_2$se)
df_2$model <- "Weighted Matched \nSample"
df_2$n <- nobs(fit1)

fit0_uw <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,data=pb2022),vocv=vcovHC)

vars <- c("Intercept","Electoral Winner")
coef <- fit0_uw[,1]
se <- fit0_uw[,2]
df_uw <- data.frame(cbind(vars,coef,se))
df_uw$cimin95 <- NA
df_uw$cimin95 <- as.numeric(df_uw$coef)-1.96*as.numeric(df_uw$se)
df_uw$cimax95 <- NA
df_uw$cimax95 <- as.numeric(df_uw$coef)+1.96*as.numeric(df_uw$se)
df_uw$cimin90 <- NA
df_uw$cimin90 <- as.numeric(df_uw$coef)-1.65*as.numeric(df_uw$se)
df_uw$cimax90 <- NA
df_uw$cimax90 <- as.numeric(df_uw$coef)+1.65*as.numeric(df_uw$se)
df_uw$model <- "Unweighted Unmatched \nSample"
df_uw$n <- nobs(fit0_uw)


fit1_uw <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner,data=pb2022mt),vcov=vcovHC)

coef_2 <- fit1_uw[,1]
se_2 <- fit1_uw[,2]
df_uw_2 <- data.frame(cbind(vars,coef=coef_2,se=se_2))
df_uw_2$cimin95 <- NA
df_uw_2$cimin95 <- as.numeric(df_uw_2$coef)-1.96*as.numeric(df_uw_2$se)
df_uw_2$cimax95 <- NA
df_uw_2$cimax95 <- as.numeric(df_uw_2$coef)+1.96*as.numeric(df_uw_2$se)
df_uw_2$cimin90 <- NA
df_uw_2$cimin90 <- as.numeric(df_uw_2$coef)-1.65*as.numeric(df_uw_2$se)
df_uw_2$cimax90 <- NA
df_uw_2$cimax90 <- as.numeric(df_uw_2$coef)+1.65*as.numeric(df_uw_2$se)
df_uw_2$model <- "Unweighted Matched \nSample"
df_uw_2$n <- nobs(fit1_uw)



df_a_w <- rbind(df,df_2)
df_a_uw <- rbind(df_uw,df_uw_2)

## Interactivee effects #####

fit2 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01, 
                    data=pb2022, weights = weight.educ.edad.pre_match), vcov=vcovHC)

vars <- c("Intercept","Low Interest","Medium Interest","High Interest")
coef <- fit2[,1]
se <- fit2[,2]
df <- data.frame(cbind(vars,coef,se))
df$cimin95 <- NA
df$cimin95 <- as.numeric(df$coef)-1.96*as.numeric(df$se)
df$cimax95 <- NA
df$cimax95 <- as.numeric(df$coef)+1.96*as.numeric(df$se)
df$cimin90 <- NA
df$cimin90 <- as.numeric(df$coef)-1.65*as.numeric(df$se)
df$cimax90 <- NA
df$cimax90 <- as.numeric(df$coef)+1.65*as.numeric(df$se)
df$model <- "Weighted Unmatched \nSample"
df$n <- nobs(fit2)

# con MatchIt

fit3 <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01,
                    data = pb2022mt, weights = weight.match.educ.edad),vcov = vcovHC)

coef <- fit3[,1]
se <- fit3[,2]
df_2 <- data.frame(cbind(vars,coef,se))
df_2$cimin95 <- NA
df_2$cimin95 <- as.numeric(df_2$coef)-1.96*as.numeric(df_2$se)
df_2$cimax95 <- NA
df_2$cimax95 <- as.numeric(df_2$coef)+1.96*as.numeric(df_2$se)
df_2$cimin90 <- NA
df_2$cimin90 <- as.numeric(df_2$coef)-1.65*as.numeric(df_2$se)
df_2$cimax90 <- NA
df_2$cimax90 <- as.numeric(df_2$coef)+1.65*as.numeric(df_2$se)
df_2$model <- "Weighted Matched \nSample"
df_2$n <- nobs(fit3)

fit2_uw <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01, 
                    data=pb2022), vcov=vcovHC)

vars <- c("Intercept","Low Interest","Medium Interest","High Interest")
coef <- fit2_uw[,1]
se <- fit2_uw[,2]
df_uw <- data.frame(cbind(vars,coef,se))
df_uw$cimin95 <- NA
df_uw$cimin95 <- as.numeric(df_uw$coef)-1.96*as.numeric(df_uw$se)
df_uw$cimax95 <- NA
df_uw$cimax95 <- as.numeric(df_uw$coef)+1.96*as.numeric(df_uw$se)
df_uw$cimin90 <- NA
df_uw$cimin90 <- as.numeric(df_uw$coef)-1.65*as.numeric(df_uw$se)
df_uw$cimax90 <- NA
df_uw$cimax90 <- as.numeric(df_uw$coef)+1.65*as.numeric(df_uw$se)
df_uw$model <- "Unweighted Unmatched \nSample"
df_uw$n <- nobs(fit2_uw)

# con MatchIt

fit3_uw <- coeftest(lm(prs_trst_w02-prs_trst_w01 ~ delta_elec_winner:pol_interesF2_w01,
                    data = pb2022mt),vcov = vcovHC)

coef <- fit3_uw[,1]
se <- fit3_uw[,2]
df_2_uw <- data.frame(cbind(vars,coef,se))
df_2_uw$cimin95 <- NA
df_2_uw$cimin95 <- as.numeric(df_2_uw$coef)-1.96*as.numeric(df_2_uw$se)
df_2_uw$cimax95 <- NA
df_2_uw$cimax95 <- as.numeric(df_2_uw$coef)+1.96*as.numeric(df_2_uw$se)
df_2_uw$cimin90 <- NA
df_2_uw$cimin90 <- as.numeric(df_2_uw$coef)-1.65*as.numeric(df_2_uw$se)
df_2_uw$cimax90 <- NA
df_2_uw$cimax90 <- as.numeric(df_2_uw$coef)+1.65*as.numeric(df_2_uw$se)
df_2_uw$model <- "Unweighted Matched \nSample"
df_2_uw$n <- nobs(fit3_uw)


df_b_w <- rbind(df,df_2)
df_b_uw <- rbind(df_uw,df_2_uw)

df_b_w$effect <- "Interactive Effects"
df_b_uw$effect <- "Interactive Effects"
df_a_w$effect <- "Main Effects"
df_a_uw$effect <- "Main Effects"

df_pb_w <- rbind(df_a_w,df_b_w)
df_pb_uw <- rbind(df_a_uw,df_b_uw)

#df_pb$n <- ifelse(df_pb$model == "Matched \nSample",nrow(pb2022mt),nrow(pb2022))


pb.plot <- ggplot(filter(df_pb_w,vars!="Intercept"),aes(model,as.numeric(coef),color=vars,group=vars)) +
  geom_hline(yintercept=0,colour="red",linetype="solid",size=.25) +
  geom_point(position = position_dodge(width = 1), size=2, aes(shape = vars)) + 
  geom_linerange(aes(ymin = cimin95, ymax = cimax95),lwd = .5, position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=cimin90,ymax=cimax90),lwd=1,position = position_dodge(width = 1)) +
  coord_flip() + theme_bw() + ylim(-0.25,0.8) +
  labs(x = "",y = "Point Estimate with 90% and 95% Confidence Intervals", caption = "", title = "",
       shape = "",color = "") + 
  facet_wrap(~factor(effect,levels=c("Main Effects","Interactive Effects"))) +
  theme(axis.text=element_text(size=12),
        axis.title=element_text(size=13),
        strip.text.y = element_text(angle=0, size=13),
        plot.caption = element_text(size=11),
        panel.grid = element_line(color = "grey70",
                                  linetype = 2)) +
  scale_color_manual(values = c("Electoral Winner" = "Black",
                                "High Interest" = "Red",
                                "Low Interest" = "Blue"),
                     breaks = c("Low Interest","High Interest")) +
  scale_shape_manual(breaks = c("Low Interest","High Interest"),
                     values = c("Electoral Winner" = 16,
                                "Low Interest" = 15,
                                "High Interest" = 17))


df_pr_w$Election <- "Presidential \nElection"
df_pb_w$Election <- "Constitutional \nPlebiscite"
df_pr_uw$Election <- "Presidential \nElection"
df_pb_uw$Election <- "Constitutional \nPlebiscite"

df_elec_w <- rbind(df_pr_w,df_pb_w)
df_elec_w$model <- paste(df_elec_w$model," \nN=",df_elec_w$n,sep="")
df_elec_uw <- rbind(df_pr_uw,df_pb_uw)
df_elec_uw$model <- paste(df_elec_uw$model," \nN=",df_elec_uw$n,sep="")

elections.plot <- ggplot(filter(df_elec_w,vars!="Intercept"),
                         aes(model,as.numeric(coef),color=vars,group=vars)) +
  geom_hline(yintercept=0,colour="red",linetype="solid",size=.25) +
  geom_point(position = position_dodge(width = 1), size=2, aes(shape = vars)) + 
  geom_linerange(aes(ymin = cimin95, ymax = cimax95),lwd = .5, position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=cimin90,ymax=cimax90),lwd=1,position = position_dodge(width = 1)) +
  coord_flip() + theme_bw() + 
  labs(x = "",y = "Point Estimate with 90% and 95% Confidence Intervals", caption = "", title = "",
       shape = "",color = "") + 
  facet_grid(factor(Election,
                    levels = c("Presidential \nElection",
                               "Constitutional \nPlebiscite"))~factor(effect,levels=c("Main Effects",
                                                                                      "Interactive Effects")),
             scales = "free_y") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),
        strip.text.y = element_text(angle=270, size=18),
        plot.caption = element_text(size=11),
        panel.grid = element_line(color = "grey70",linetype = 2),
        strip.text = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.position = "bottom") +
  scale_color_manual(values = c("Electoral Winner" = "Black",
                                "High Interest" = "Red",
                                "Low Interest" = "Blue"),
                     breaks = c("Low Interest","High Interest")) +
  scale_shape_manual(breaks = c("Low Interest","High Interest"),
                     values = c("Electoral Winner" = 16,
                                "Low Interest" = 15,
                                "High Interest" = 17))


elections.plot_uw <- ggplot(filter(df_elec_uw,vars!="Intercept"),
                         aes(model,as.numeric(coef),color=vars,group=vars)) +
  geom_hline(yintercept=0,colour="red",linetype="solid",size=.25) +
  geom_point(position = position_dodge(width = 1), size=2, aes(shape = vars)) + 
  geom_linerange(aes(ymin = cimin95, ymax = cimax95),lwd = .5, position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=cimin90,ymax=cimax90),lwd=1,position = position_dodge(width = 1)) +
  coord_flip() + theme_bw() + ylim(-0.18,0.8) +
  labs(y = "",x="") + 
  #labs(x = "",y = "Point Estimate with 90% and 95% Confidence Intervals", caption = "", title = "",
  #     shape = "",color = "") + 
  facet_grid(factor(Election,
                    levels = c("Presidential \nElection",
                               "Constitutional \nPlebiscite"))~factor(effect,levels=c("Main Effects",
                                                                                      "Interactive Effects")),
             scales = "free_y") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),
        strip.text.y = element_text(angle=270, size=18),
        plot.caption = element_text(size=11),
        panel.grid = element_line(color = "grey70",linetype = 2),
        strip.text = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.position = "bottom") +
  scale_color_manual(values = c("Electoral Winner" = "Black",
                                "High Interest" = "Red",
                                "Low Interest" = "Blue"),
                     breaks = c("Low Interest","High Interest")) +
  scale_shape_manual(breaks = c("Low Interest","High Interest"),
                     values = c("Electoral Winner" = 16,
                                "Low Interest" = 15,
                                "High Interest" = 17))


df_elec_uw$weight <- "Unweighted"
df_elec_uw$model <- gsub("Unweighted","",df_elec_uw$model)

df_elec_w$weight <- "Weighted"
df_elec_w$model <- gsub("Weighted","",df_elec_w$model)

df_elec_full <- rbind(df_elec_uw,df_elec_w)


ggplot(filter(df_elec_full,vars!="Intercept"),
                         aes(model,as.numeric(coef),
                             color=factor(vars,levels = c("Electoral Winner","Low Interest","Medium Interest",
                                                          "High Interest")),
                             group=factor(vars,levels = c("Electoral Winner","Low Interest","Medium Interest",
                                                          "High Interest")))) +
  geom_hline(yintercept=0,colour="red",linetype="solid",size=.25) +
  geom_point(position = position_dodge(width = 1), size=2, aes(shape = vars)) + 
  geom_linerange(aes(ymin = cimin95, ymax = cimax95),lwd = .5, 
                 position = position_dodge(width = 1)) +
  geom_linerange(aes(ymin=cimin90,ymax=cimax90),lwd=1,
                 position = position_dodge(width = 1)) +
  coord_flip() + theme_bw() + 
  labs(x = "",y = "Point Estimate with 90% and 95% Confidence Intervals", caption = "", title = "",
       shape = "",color = "") + 
  facet_grid(factor(Election,
                    levels = c("Presidential \nElection",
                               "Constitutional \nPlebiscite"))~factor(effect,
                                                                      levels=c("Main Effects",
                                                                               "Interactive Effects")) + 
               weight,scales = "free_y") +
  theme(axis.text=element_text(size=18),
        axis.title=element_text(size=18),
        strip.text.y = element_text(angle=270, size=18),
        plot.caption = element_text(size=11),
        panel.grid = element_line(color = "grey70",linetype = 2),
        strip.text = element_text(size = 18),
        legend.text = element_text(size = 18),
        legend.position = "bottom") +
  scale_color_manual(values = c("Electoral Winner" = "Black",
                                "High Interest" = "#F8766D",
                                "Medium Interest" = "#00BA38",
                                "Low Interest" = "#619CFF"),
                     breaks = c("Low Interest",
                                "Medium Interest",
                                "High Interest")) +
  scale_shape_manual(breaks = c("Low Interest",
                                "Medium Interest","High Interest"),
                     values = c("Electoral Winner" = 16,
                                "Low Interest" = 15,
                                "Medium Interest" = 18,
                                "High Interest" = 17))

```
